Renormalization group contraction of tensor networks in three dimensions 
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We present a new strategy for contracting tensor networks in arbitrary geometries. This method 
is designed to follow as strictly as possible the renormalization group philosophy, by first contracting 
tensors in an exact way and, then, performing a controlled truncation of the resulting tensor. We 
benchmark this approximation procedure in two dimensions against an exact contraction. We then 
apply the same idea to a three dimensional system. The underlying rational for emphasizing the 
exact coarse graining renormalization group step prior to truncation is related to monogamy of 
entanglement. 
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Computing analytically the properties of a quantum 
model is, in general, not possible. It is then necessary 
to resort to classical simulations that rely on approxi- 
mation techniques. From a quantum information per- 
spective, the presence of a small amount of entanglement 
in the system has been identified as the key ingredient 
that allows for efficient classical simulations. This has 
been exploited in one dimension (ID) in a series of al- 
gorithms based on a description of the system as a ten- 
sor network. Following this approach, condensed matter 
systems can be simulated using Matrix Product States 
(MPS, as the testbed for the variational procedure 

of the density-matrix renormalization group (DMRG, Q) 
to study ground states of local ID Hamiltonians. Beyond 
the ID setting, the computation of physical magnitudes 
using tensor networks is limited by the numerical effort 
necessary to perform the contraction of the tensor net- 
work, i.e. to sum all its indices. The efficiency for per- 
forming this task is limited by the area law scaling of the 
entanglement entropy in the system To overcome 

this problem, several strategies aim at finding the best 
possible approximation to the contraction of tensor net- 
works after identifying the relevant degrees of freedom of 
the system 

Let us briefly recall the key elements of the tensor net- 
work representation. Given a quantum state of N par- 
ticles = ^ c'l''"'*" |ii . . . «Ar) its coefficients can be 
represented as a contraction of local tensors c*i' - '*« = 
tr(A^'*i . . . A^'*"), where each local tensor A' carries a 
physical index ij and ancillary indices (which are not 
written) that get contracted according to a prescribed 
geometry. The rank of these ancillary indices, that we 
shall call Xj controls the amount of entanglement which 
is captured by the tensor representation. If the tensors 
are simple matrices on a line, the tensor network is called 
MPS [lill]. Other possible geometries are regular squared 



grids in any dimensions that correspond to PEPS [ll|, 
and tree-like structures that go under the name of TTN 
[H and MERA [13]. 

In this letter we propose a new strategy to contract 
tensor networks in general geometries, that we shall il- 



lustrate in detail for PEPS in 2D and 3D. The method is 
based on following as strictly as possible the renormaliza- 
tion group (RG) philosophy. First, an exact contraction 
of a set of local tensors is performed that produces a 
coarse grained tensor of larger rank. Subsequent con- 
tractions would made the rank of the effective tensors 
to scale following a law dictated by the geometry of the 
tensor network. For instance, in the case of PEPS the 
rank of the coarse grained tensors would grow following 
an area law. It is then necessary to perform a truncation 
that faithfully retains only relevant degrees of freedom. 
To achieve this truncation, our method makes a series of 
Schmidt decompositions that cast the relevant informa- 
tion onto a renormalized tensor. 

The strategy we present here differs from a previous 
renormalization group inspired proposal 14|-|l6l|. There, 
the original tensor is first truncated and then contracted 
efficiently. Instead, we first contract exactly tensors at a 
larger numerical cost, and then truncate. This procedure 
can be made exact in ID I7|. In more dimensions, this 
idea entails a trade off between numerical computation 
speed and precision. Our proposal relies on the idea that 
in higher dimensions, e.g. 3D, monogamy of entangle- 
ment makes every degree of freedom to have a reduced 
amount of entanglement with each neighbor. Long dis- 
tance correlations emerge from the multiplicity of possi- 
ble paths connecting local degrees of freedom. Therefore, 
a tensor network with small rank x is already a good ap- 
proximation to a 3D system. The fact that good tensor 
network representations only need a small x makes viable 
our proposal for exact contractions followed by trunca- 
tion. 

Let us recall that, in 2D settings, some variants of con- 
traction schemes for PEPS have been analyzed [HI, 

To our knowledge, PEPS renormalization methods 
have not been applied to 3D systems, which nevertheless 
have been studied using cluster states [l8| or string bond 
states [l^. For 3D classical systems other renormaliza- 
tion algorithm have been proposed 2ol 21|. 

RG contraction in 2D. — Let us illustrate our method 
for contraction of tensor networks with the computa- 
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tion of the norm of a state in 2D. We start by con- 
sidering the norm of the original state that is the 
folded lattice made of tensors of the form = 

^^^*(m,n)j ^ ^j^^^.^, (-^^ ^) Yshcis a position in 

the 2D lattice and physical indices i at this point have 
been summed over. Each tensor £'('">") has four fur- 
ther ancillary indices of rank [not written explicitly) 
pointing to its 2D neighbors. The global contraction of 
the 2D tensor network corresponds to the computation 
of the norm (^/^IV^) = Tr (i^^^^^) . . . £;(^'^)) , where Tr is a 
generalized trace that contracts all ancillary indices. 

The stages in our contraction strategy go schemat- 
ically as follows. We start with the tensor 
which carries 4 indices of rank x^- Wc first contract 
four adjacent tensors to build a plaqucttc; = 

j^{rn,n) j^{m,n+l) j^(m+l,n) j^(m+l,n+l) ^ Now P'™'"' Carries 

a total of 8 open indices of rank x^, showing the area law 
increase of indices naturally associated to the contrac- 
tion of tensors in a 2D square lattice. This plaqucttc 
tensor will be approximated by a renormalized tensor on 
a coarse grained lattice S^™'"' p(™.") ^ where we shall 
truncate to a tensor of 4 indices of rank . The global 
contraction symbolically reads 

(VI V) = TV {{E}) = Tr ({P}) « Tt ({£}) , (1) 

which can be graphically represented as 




where tensors E arc represented by nodes connected to 
its neighboring tensors with lines that represent ancil- 
lary indices. After a renormalization step, we recover 
the same lattice structure with renormalized tensors. 

Let us now discuss in more detail each step in the 
renormalization group contraction we have just sketched. 
Though computationally expensive, the initial exact con- 
traction {E} — >• {P} has two main advantages. First, an 
exact contraction of the tensors in a plaquette retains all 
degrees of freedom, making this method specially suitable 
for systems with frustration at the plaquette scale. Fur- 
thermore, by an appropriate renormalization of the vir- 
tual bonds of the P tensors we recover the original lattice 
structure, modulo an increase in the rank of ancillary in- 
dices. When completing the strategy with a truncation, 
it will then be possible to perform a systematic iteration 
of the same procedure, resulting in the global contraction 
of the tensor network. The method is thus a way to scape 
from the area law restriction at the cost of a truncation 
of the renormalized tensors. 



Once we have constructed the renormalized tensor 
, we need to perform a controlled truncation. Let 
us modify the notation to make this step clearer, by drop- 
ping the site label and introducing explicitly the ancil- 



lary indices P^^ 



1) 



Each ancillary index a 



is the result of combining the two ancillary indices that 

were attached to tensors E pointing towards an adjacent 
plaquette. To be concrete, we call ai the index connect- 
ing the plaquette to its upper neighbor. This index is 
the combination of two indices from the original tensors 
a\ = (3i. Thus, each a is a combined index of rank 
X*- We then select the index a\ and separate it from the 
rest of indices of the plaquette by means of a singular 
value decomposition 



(2) 



where U and V are unitary matrices and arc the eigen- 
values in the decomposition. It is convenient to include 
the squared root of the eigenvalues into the unitary ma- 



trix C/ai/n as follows 



(3) 



We repeat this process for each index in order to obtain 
the four matrices W\, W2, W3 and W4 as follows 
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After these four successive decompositions, we can con- 
struct a new tensor P = W{~^W2^W^^W^^P, where we 
have used pseudo inverse matrices. 

We proceed now to perform a renormalization of the 
plaquette P using the set of tensors Ws. The key obser- 
vation here is that, along each direction, the plaquette is 
surrounded by other plaquettes decomposed in a similar 
way. We assume that a plaquette P = W1W2WSW4P has 
neighboring unitaries W^^W^^W^ and W4. The renor- 
malization consists on the truncation of the degrees of 
freedom between plaquettes through the unitaries W and 
W . We first take matrices Wi and Wg to form tensor 
Ri,i = WgW^i which can be decomposed using a SVD into 
Ps i = tlYiV . The tnmcation corresponds to only keep- 
ing the largest x^ eigenvalues from matrix S to generate 
the new decomposition i?3^i Ri /3/1, where each matrix / 
is the projection to W onto the relevant eigenvalue sub- 
space. We repeat this procedure for each neighboring 
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pair W and W . We then apply the four matrices Ik to 
P to obtain our final truncation E = /1/2/3/4P 



o — o 



Q 




This truncation strategy that delivers the renormalized 
tensor E can be represented as 

P = Tr{EEEE) = WiW2W3WiP « hhhhP = E 

and reduces the size of the plaquette tensor to for each 
bond. 

The procedure described above is computationally de- 
manding compared to variational methods or other renor- 
malization techniques. In order to make the singular 
value decomposition in Eq. [5] it is necessary to make 
the contraction of a plaquette with its own adjoint 



(4) 

to obtain each of the Wk ■ This strategy outperforms the 
naive use of singular value decomposition algorithms. 

As mentioned above, previous proposals for the renor- 
malization of tensors networks rely on a decomposition of 
the local tensors before producing a truncation 14-16]. 
Such a decomposition allows to achieve a higher bond di- 
mension after the original truncation. Other approaches 
work with square plaquettes that include the physical in- 
dices [22]. It is also important to notice that the method 
proposed here is not a variational procedure, though it 
remains numerically stable. 

Validation in 2D.- We validate our approach showing 
results for the renormalization of a tensor network rep- 
resenting the ground state of the Ising Hamiltonian with 
a transverse magnetic field in a square lattice 



H 



E 



(5) 



where {i,j} are neighboring sites on finite or infinite lat- 
tices. This model displays a phase transition at a critical 
field he ~ 3.04 for the infinite square lattice. 




FIG. 1. For an Ising model with transverse field in a fi- 
nite 2D square lattice of size 6 x 6, we plot the error of 
the magnetization between the exact contraction and our 
renormalization group contraction of the same PEPS, that 
is mz(x, exact) — mz{x, RG), for difi^erent values of x- Inset: 
Same for a 12 x 12 lattice with x = 2. 



The preparation of the tensor network can be per- 
formed using an imaginary time evolution 11, l^], IV'/) = 



e~ IV'o)- We first proceed to use a Trotter approxima- 
tion to this Euclidean evolution. At each step we have to 
apply the evolution operator T to a pair of neighboring 

and truncate the 



tensors 6„'^ 



and [Al^ 



to the lattice bond 



new tensors [^^^^ 
dimension. 

In order to validate our strategy, we show in Fig[T] the 
error of the magnetization ruz in a finite square lattice of 
size 6x6 and 12 x 12, with our RG contraction against 
the exact contraction of the very same tensors, that is 
|mz(x, exact) — mz{Xi RG)\. The error introduced by 
the RG contraction is of the order of 10^"^ at the crit- 
ical point. Let us note that this error is smaller for x = 2 
since the RG casts x^^ numbers to x*. It is thus expected 
that the faithfulness of RG is better for small x- Obvi- 
ously, the X = 4 tensor is in itself a better representation 
of the state. 

3-D quantum systems.- The contraction strategy for 
tensor networks we have presented in 2D can be extended 
for a 3D square lattice. The renormalization is performed 
over plaquettes of eight tensors forming a cube. We use 
a singular value decomposition to decouple at each step 
four ancillary indices, and obtain the unitaries corre- 
sponding to each orthogonal direction. This set of six 
matrices Wk is combined with unitaries W^. associated 
to neighboring plaquettes, to produce six truncated ma- 
trices Ik- We perform a similar operation along each 
direction to recover the renormalized tensors with bond 
dimension x^- To overcome limitations in computational 
resources, we use an operation similar to the trick in |3] 
and obtain each of the unitaries Wi. For a single direc- 
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FIG. 2. Magnetization vs. local field h for the Ising 
model with transverse field in an infinite 3D square lattice. 
The critical point is found at he ~ 5.29. The rank of the 
tensor network is x = 2. Inset: Magnetization rUx vs. h for 
different values of Xeff- 




FIG. 3. Renormalization of hexagonal and triangular lattices, 
forming plaquettes with the tensors inside the shaded areas. 



be contracted as in Fig|3]so as to recover a rescaled ver- 
sion of the original lattice. This opens the possibility of 
studying frustrated models in 3D using PEPS technology. 
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tion, this renormalization is represented as 




(6) 



It is possible to prepare the original state with a tensor 
network using some specific x- Yet, after a first contrac- 
tion is made, we can retain renormalized tensors with a 
larger effective Xeff- We have checked that for x = 2, 
stable results are obtained for Xeff ~ 5. 

We can now apply our approach to 3D systems to com- 
pute the magnetization of Ising model Eql5]for an infinite 
3D square lattice, by iterating the above procedure. The 
state is prepared using a minimal environment to stabi- 
lize the Trotter evolution. The results for the magneti- 
zations ruz and for x — ^re presented in Figl2] A 
quantum phase transition is detected at a critical point 
located around he ~ 5.29 (series expansions detect a crit- 
ical point at he = 5.14 24|). We have also computed 



the equivalent phase transition in 2D, using again only 
X = 2 tensors. In the 2D case, our KG contraction finds 
a transition at he ~ 3.25 to be compared with the re- 
sult he = 3.04 coming from other more precise methods. 
This hints at the fact that RG contraction is a better 
approximation in 3D than in 2D. 

Conclusions. — We have presented a novel, RG in- 
spired strategy to contract tensor networks that delivers 
good results in a 3D simulation of the phase transition in 
the quantum Ising model. Our scheme can be extended 
to hexagonal and triangular lattices, where a renormal- 
ization step can be performed choosing the plaquettes to 
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